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We introduce the numerical linked cluster (NLC) expansion as a controlled numerical tool for 
the study of the many-body localization (MBL) transition in a disordered system with continuous 
non-perturbative disorder. Our approach works directly in the thermodynamic limit, in any spatial 
dimension, and does not rely on any finite size scaling procedure. We study the onset of many-body 
delocalization through the breakdown of area-law entanglement in a generic many-body eigenstate. 

By looking for initial signs of an instability of the localized phase, we obtain a value for the critical 
disorder, which we believe should be a lower bound for the true value, that is higher than current 
best estimates from finite size studies. This implies that most current methods tend to overestimate 
the extent of the localized phase due to finite size effects making the localized phase appear stable 
at small length scales. We also study the mobility edge in these systems as a function of energy 
density, and find that our conclusion is the same at all examined energies. 

PACS numbers: 


Introduction — The eigenstate thermalization hypoth¬ 
esis (ETH) is a powerful statement relating observables 
of the high energy eigenstates of a quantum many-body 
system with their thermal expectation values UM- How¬ 
ever, this principle can be violated in certain systems 
with strong enough disorder, where even the l^h energy 
eigenstates possess only local entanglement UI3|- An¬ 
derson localization is a one-body example of this. An 
area of key interest is how far this localization persists 
in a many-body state in presence of interactions Q. At 
what point are interactions strong enough that the local¬ 
ization is destroyed and the system obeys ETH? This 
is the problem of many-body localization (MBL) transi¬ 
tion, which is a topic of active research both theoretically 
and experimentally dill- 


The surge of interest in many-body localized systems 
has motivated many numerical studies. Most studies 
have focused on exact diagonalization or Lanczos meth¬ 
ods which are able to address both sides of the tran¬ 
sition in small systems (siUisjl . However, since much 
about this phase transition is still not well understood, 
extension of finite size results to the thermodynamic limit 
can prove difficult. We would like to examine this phase 
transition using expansion methods, which provide an 
alternate way of addressing the thermodynamic limit. 
While standard perturbative series expansions are very 
p ower ful (idl - liil , they suffer from small energy denomi¬ 
nators in models with continuous non-perturbative disor¬ 
der. Thus, we turn to the numerical linked cluster (NLC) 
expansion Ublj, which does not suffer from this prob¬ 
lem of small energy denominators. 

In this paper, we provide evidence that for a prototypi¬ 
cal model of MBL, approaching the critical disorder from 
the localized side, the localized phase actually becomes 
unstable only at increasingly long length scales inacces¬ 
sible to most numerical techniques. This implies that 
finite size numerical studies on the MBL transition tend 


to overestimate the extent of the localized phase. 

Model — The system we study explicitly is the spin-1/2 
Heisenberg Hamiltonian with random fields along the z 
direction, 


n = J2h.sf + J2s^■s„ ( 1 ) 


on the Id chain, where the sum {i,j) is over adjacent 
pairs. The random field is picked from a uniform distri¬ 
bution hi G [—h, h]. This is one of the simplest models to 
study many-body localization ^ and has been studied 
numerically in great detail |32H35l l52ll53l| . At low h, the 
system is in a thermalizing phase obeying the ETH, while 
at high h, the system is in the localized (MBL) phase. 

To identify these different phases, we focus on the en¬ 
tanglement properties of the eigenstates. The typical 
measure for entanglement in a pure state bipartitioned 
into two parts A and B is the von Neumann entropy, 
defined for some state |d>) as 

s(l^)) =-Tr(pAlnpA,) (2) 

where pA = Tr^ jT*) ('Ll is the reduced density matrix, 
obtained by tracing over all external degrees of freedom 
from the density matrix. 

A typical eigenstate in an ETH obeying system will 
exhibit thermal volume law entanglement. The entangle¬ 
ment entropy will approach the classical thermal entropy 
(required by ETH) and scale with the volume of the re¬ 
gions A and B. In the localized phase, the eigenstates 
will instead obey an area-law, scaling with the area be¬ 
tween A and B. This can be understood by regarding 
them as simultaneous eigenstates of many local opera¬ 
tors Q: only due to mixing contained in operators near 
the boundary will one get contributions to the entangle¬ 
ment, which therefore grows with the area of the bipar¬ 
titioning. 
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Numerical Linked Cluster Expansions — NLC is simi¬ 
lar to perturbative series expansions in that interactions 
within clusters of increasing sizes must be considered, 
but rather than perturbatively treating each interaction 
within a cluster, we solve them numerically, typically by 
exact diagonalization. Our treatment of disorder in NLC 
is different from usual ^ M\, allowing us to deal with 
continuous non-perturbative disorder. The procedure is 
outlined briefly here. 

Let N be the order to which we wish to do the calcula¬ 
tion. The order of the calculation defined as the number 
of spins in the largest cluster considered. We identify a 
finite size region of the infinite system to work with. For 
the chain, this is simply a 2{N — 1) length chain, with a 
bipartitioning cut in the middle. Each of the sites i is as¬ 
signed a field hi, which is held fixed until the calculation 
is complete. The reason for this choice of system size is 
so that the results remain correct for the infinite system, 
to the desired order, as explained later. 

We define a cluster c to be a set of sites. A Hamiltonian 
T-lc can be obtained considering only the spins in c, and 
diagonalized numerically to obtain the eigenstates {|a'^)} 
with eigenvalues {e° } labeled by a. Our quantity of in¬ 
terest, the eigenstate averaged entanglement entropy, can 
then be calculated as 


S(.c)=J2^s{\a^)), (3) 

CK 

where Z = is the normalization factor and 

/3 = 1 /T is the inverse temperature. 

The entropy for the infmite lattice C can be expressed 
as a sum over the weight S{c) of all clusters c that can be 
embedded in the lattice: S{jC,) = The weight 

of a cluster is then defined recursively by the principle of 
inclusion and exclusion: 


S{c) = S{c) - ^ S{c’). (4) 

c'Cc 


One can show that only connected clusters which cross 
the boundary can have a nonzero weight. First, if a 
cluster does not cross the boundary there can obvi¬ 
ously be no entanglement in it or its subclusters, so the 
weight is trivially zero. Second, proving that only con¬ 
nected clusters can contribute simply amounts to prov¬ 
ing that S obeys the linked cluster property, that is, for 
a cluster with two disconnected components ci and C 2 , 
S'(ci U C 2 ) = S'(ci) + S'(c 2 ). This follows from the fact 
that 'HciUc 2 = hici (B'Hc 2 - Thus, we must simply consider 
connected clusters of up to size N that have sites on both 
sides of the partition. Our finite size representation was 
chosen to contain all the necessary clusters of the infinite 
system up to order N. The count for the clusters cross¬ 
ing the boundary scales with the area thus guaranteeing 
an area-law as long as NLC converges. 


Using this, we can obtain a series a„ whose sum gives 
the total eigenstate averaged entanglement entropy per 
unit area ^area, 


e — 


N 

■E' 

n=0 


= E 

c,\c\—n 


where Lcut is simply 1 for the chain. Finally, the entire 
calculation must be repeated for different realizations of 
{hi} to obtain a disorder averaged value for a„ 
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The NLC scheme used here is slightly different from 
what has typically been used for random systems in the 
past, where different embeddings of the same graph are 
treated identically and one does not need a consistent 
finite system [3 t 5^- The more standard scheme has 
been _Mplied to study MBL systems with discrete disor¬ 
der m, which allows one to perform disorder averaging 
before subgraph subtraction. When there is continuous 
disorder, partial disorder averaging over a finite number 
of realizations means that the linked cluster property is 
only approximate, and thus large errors will build up at 
high orders. In our approach, the linked cluster prop¬ 
erty is guaranteed and one is free to average over many 
realizations of the system. This is a key aspect of our 
calculation which allows us to treat disordered systems 
with continuous non-perturbative randomness. 

Does it converge? — We first examine T = 00 . If en¬ 
tanglement satisfies a thermal volume-law, interpreting 
n as a proxy length scale [s^, we expect a„ to eventu¬ 
ally saturate to the volume-law constant ln(2)/2 at high 
enough n [^. We should note that our model (Eq[T]) 
does not possess a strongly thermalizing regime, due to 
the integrability at h = 0, and therefore we do not yet 
see this saturation to the thermal value within our range 
of n [ 5 ^. In the localized phase, the additional entangle¬ 
ment due to the addition of one site far away from the 
cut should become exponentially small with distance, so 
we expect a„ to decay exponentially to 0 once n is larger 
than some localization length f. We define the MBL 
phase in our study to be one in which the sum of On 
converges exponentially. 

Fig[T]shows On for a range of h values. To estimate con¬ 
vergence or divergence, a linear extrapolation to 1/n = 0 
can be performed. If the extrapolation predicts Ooo > 0, 
we argue that this corresponds to a breakdown of area- 
law. Although we expect a„ to eventually go to 0 expo¬ 
nentially in the localized phase, this would only happen 
when our cluster sizes are much greater than the local¬ 
ization length scale. This can be more clearly shown by 
examining the ratio of terms, which we will discuss next. 
Note that the case of an area-law with logarithmic cor¬ 
rections would correspond to one where a„ heads linearly 
to 0, which we consider in this analysis the boundary be¬ 
tween convergence and divergence. 

Let us define the ratio of the (disorder averaged) se¬ 
ries terms = a„/a„_i. In the MBL phase, we can say 
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FIG. 1: The nth order area-law contribution a„ (Eq.[5} at T = 
oo. Near he, data has been averaged over more than 3 x 10® 
disorder realizations of the chain, and error bars show the 
standard error of the mean. Dashed lines are a demonstration 
of the linear extrapolation to a^o by fitting the last 4 terms 
in the series. 



1/n 


FIG. 2: Plot of the ratios r„ = On/un-i versus 1/n. Also 
shown is the line exp( —1/n), above which we argue the ex¬ 
pansion cannot converge exponentially. 


more about the overall trend of r„. Again interpreting 
n as a proxy length scale, for large n, we expect an to 
decrease exponentially with potentially power-law pref¬ 
actors. The leading contributions at n ^ should be of 
the form a„ = Cn~^ exp{—n/^), where k is some posi¬ 
tive number, ^ is the localization length, and C is some 
arbitrary constant Q. Therefore, in the large n limit, 
discarding terms smaller than 1/n, we expect 

rn = an/on-i = (1 - A:/n) exp(-l/^). (6) 

Therefore, plotting r„ versus 1/n, r„ should approach 
T'exi — exp(—1/^) from below, with a slope of —k. How¬ 
ever, near the transition ^ can become very large and we 
do not actually see this behavior within our range of at¬ 
tainable n. We can, however, predict whether this kind 
of exponential convergence is possible given the behavior 
of the series at finite n. 


Fig shows the behavior of for our range of h and 
n. The trend seems to be for to increase steadily 
with n (although eventually r„ must approach 1 in the 
delocalized phase). If the system is in the MBL phase and 
the series is to converge exponentially, we expect that at 
some n a„ will begin decreasing exponentially and 
Tn will begin heading towards Too with slope —k. Barring 
bizarre behavior such as increasing to some high value 
and then suddenly decreasing before finally increasing 
again towards Too, this places a restriction on what 
can be when a„ begins its exponential decay. Because 
the slope of the approach is negative or 0, < Too- 

But also n ^ which means this decay can only begin 
occurring when 


Tn <roc = exp(-l/^) < exp(-l/n). 


( 7 ) 


Therefore, once has increased above exp(—1/n), a„ 
cannot converge exponentially. This is not a rigorous 
claim, but should be valid as long as a„ behaves in a 
regular manner. This clearly shows (in Fig j^]) that the 
series for h = 4.0 cannot converge exponentially, and thus 
is not in the MBL phase. However, the series for h > 4.5 
are still within this region, and thus may diverge or con¬ 
verge. Hence, our result should serve as a lower bound, 
with our best estimate being at he = 4.5 ±0.1. Going to 
higher order in NLC can further refine this value. 

Discussion — A way of viewing our result is that 
we are seeing an instability to thermalization of an 
almost-localized regime [^ . That is (in Fig[T]), initially 
On acts quite localized in that it is much smaller than the 
thermal value and is getting smaller as n is increased, but 
may start increasing at higher n, signaling the “onset” of 
thermalization. This onset of thermalization moves to 
higher n as h is increased, but goes beyond our range of 
accessible n after h = 3.5. However, by looking in a sen¬ 
sitive way for initial signs of an instability, we are able 
to place a lower bound for he at 4.5 ± 0.1. 

To understand how our analysis is more sensitive to 
this transition than other methods, let us focus on the 
entanglement per unit volume S'voi- •S'voi decreasing with 
system size is often associated with area-law entangle¬ 
ment and therefore localization [^. In our study, S'voi 
for a system of size N would correspond to the quan¬ 
tity Svoi = {^/N)J2n=o^^- already 

“turned up” and was increasing (clearly thermalizing), 
Svoi would not begin increasing until a„ had increased 
above the mean of all the previous terms in the series. 
Our analysis predicts this upturn, which itself would pre¬ 
cede estimates from finite size systems using Svoi- 

Other methods, which are more focused on seeing the 
full onset of thermalization, do not observe the transition 
near our bound 33- 35|. Results from Lanezos on systems 
of up to 22 sites with finite size scaling show evidence for a 
transition at /ic ~ 3.7 33|. However, the scaling exponent 
^ ~ (/i — hc)~'^ obtained from finite size scaling strongly 
violates the Harris-Chayes bound in Id [13, evidence 
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FIG. 3: The phase diagram with (scaled) energy and disor¬ 
der on the axes, identified by interpolation of the intercept 
method (Fig [T|) up to order 10. Error bars represent confi¬ 
dence in our interpolation. Also shown is the he obtained by 
finite size scaling of Lanezos results from Ref (error bars 
not shown). 


that perhaps they are still far from the true critical point. 
This implies that finite size effects are significant, even in 
the systems accessible to Lanezos, and some corrections 
to the finite size scaling is needed. These effects cause 
an overestimate of the stability of the MBL phase, which 
actually becomes unstable to thermalization earlier only 
at much longer length scales. Note that an interesting 
alternate possibility is the existence of an intermediate 
phase between the ETH and MBL phase, with neither 
thermal nor area-law entanglement 63“66|, which we do 
not pursue further. 

This onset of thermalization at high order is what 
one would expect from a long lengthscale delocalization 
mechanism. In studies of this transition by a renormal¬ 
ization group approach, one also finds that the transi¬ 
tion is driven by rare metallic inclusions [b^, [b^. Near 
the transition, these are rare enough that small systems 
look localized, but actually become thermalizing at long 
length scales. 

The mohility edge — Finally, we can observe the tran¬ 
sition at different energy windows by varying /3 = 1/T 
in Eq. [3] of our NLC calculation, thus probing states at 
a given energy defined by the thermal ensemble. Fol¬ 
lowing the same arguments as at T = oo, we can ob¬ 
tain estimates for he at a given temperature or en¬ 
ergy density. Fig [3] shows our estimates for various (3 
values, with the scaled energy e on the vertical axis: 
e = (F F'min)/(F/niax .^min); where T/eain and F^ax ‘^'I’e 
the lowest and highest energies in the energy spectrum. 
The shape of our estimates are very similar to those ob¬ 
tained in previous numerical calculations 32, |^, along 
with the slight asymmetry expected around e = 1/2 [z3- 
As with the case at T = oo, we find our estimates are 
consistently higher than previous numerical calculations. 

There is much debate on whether a mobility edge exists 


in the thermodynamic limit. Numerical results suggest 
the presence of such an ed^ 32, 33, but there are 
also arguments against it [b7|. While our phase diagram 
shows a similar shape as previous numerical studies, our 
analysis gives a lower bound for he, and hence does not 
negate the claims of the absence of a mobility edge. If the 
transition actually occurs at a single he for all energies, 
the fact that our estimates are lower away from the center 
of the spectrum would indicate that much larger length 
scales would be needed to observe delocalization and that 
finite size effects would be much stronger in those regions. 

Conclusions — In conclusion, we have studied the MBL 
transition in the random field Heisenberg model using 
NLC expansions. We focus on the breakdown of the area- 
law of entanglement in the eigenstates of the Hamilto¬ 
nian. Our approach works directly in the thermodynamic 
limit and does not rely on any finite size scaling. By look¬ 
ing for signs of instability in the MBL phase, we are able 
to estimate a lower bound for the critical disorder in the 
thermodynamic limit. At all energies examined, our he 
estimates are consistently higher than found by finite size 
studies. This implies that numerical methods which look 
for the full onset of thermalization tend to overestimate 
the extent of the MBL phase, which actually becomes 
unstable earlier but only at much longer length scales. 
Near the transition, finite size effects are significant and 
hence caution must be taken when relating to the infinite 
system. 

Our result can be readily verified by cold atoms ex¬ 
periments, which are able to present very well character¬ 
ized systems [isl - l^ . If the Id random field Heisenberg 
model is experimentally realized, measurements of the 
critical disorder for large systems should lie above our 
estimate. We may also extend our result for other simi¬ 
lar models of the disorder driven MBL transition, where 
delocalization also occurs over a long lengthscale [68], l69j , 
and suggest that the true critical point would be higher 
than finite size scaling estimates from small systems. 
Also of interest are quantum chaotic Wannier-Stark sys¬ 
tems (nUii, which are experimentally accessible and 
possess a localization-delocalization transition in the ab¬ 
sence of disorder. 

We would like to thank David Huse for many valuable 
discussions. This work is supported in part by NSF grant 
number DMR-1306048. 
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